######## MAIN MEETING-LEVEL ANALYSIS ########

source("./code/loadPackages.R") # Install and load all necessary packages

## Load functions to make coefficient plots
source("./code/makeCoefPlots.R")

## Load meeting data
nscd_all = fread("./data/meetingData.csv")
nscd_all$admin = factor(nscd_all$admin, levels=c("Truman", "Eisenhower", "Kennedy", "Johnson", "Nixon", "Ford", "Carter", "Reagan"))
nscd_all$date = as.Date(nscd_all$date)


## Set up variables
outcome = "nConfAdv"
outDiff = "nConfAdv-nCoopAdv"
controls = " + nAttend + numDef + numInt + numMil + numState + 
              log(totDip+1) + log(totInt+1) + log(totMil+1) +  
              log(usmidschallenge5+1) + cinc + formal +
              mention.Economy + mention.Europe + mention.Asia + 
              mention.Defense + mention.China + mention.Other + 
              mention.Intelligence + mention.Middle.East + mention.Strategic.Forces + 
              mention.Organization + mention.USSR + mention.Policy + 
              mention.Americas + mention.International.Institutions + mention.Diplomacy + 
              mention.Vietnam + mention.Africa + mention.Arms.Control + 
              mention.North.Africa + mention.Latin.America + factor(admin)"
controlsNoAdmin = gsub("\\+ factor\\(admin\\)", "", controls)


## Emergence model
meanFit1a = glm(paste(outcome, paste0("~ meanHawk", "+ formal + factor(admin)"), sep=" "),  
                data=nscd_all, family="poisson")
summary(meanFit1a)

meanFit1b = glm(paste(outcome, paste0("~ meanHawk", "+", controls)),  
                data=nscd_all, family="poisson")
summary(meanFit1b)

meanFit1c = lm(paste(outDiff, paste0("~ meanHawk", "+ formal + factor(admin)"), sep=" "),  
               data=nscd_all)
summary(meanFit1c)

meanFit1d = lm(paste(outDiff, paste0("~ meanHawk", "+", controls)),  
               data=nscd_all)
summary(meanFit1d)


## Leader model
leaderFit1a = glm(paste(outcome, paste0("~ presHawk", "+ formal"), sep=" "),  
                  data=nscd_all, family="poisson")
summary(leaderFit1a)

leaderFit1b = glm(paste(outcome, paste0("~ presHawk", "+", controls, "-factor(admin)", sep=" ")),  
                  data=nscd_all, family="poisson")
summary(leaderFit1b)

leaderFit1c = lm(paste(outDiff, paste0("~ presHawk", "+ formal"), sep=" "),  
                 data=nscd_all)
summary(leaderFit1c)

leaderFit1d = lm(paste(outcome, paste0("~ presHawk", "+", controls, "-factor(admin)", sep=" ")),  
                 data=nscd_all)
summary(leaderFit1d)


# Likelihood ratio test for group mean model to test for emergence
meanFit1b_null = glm(paste(outcome, "~ ", controls, sep=" "),  
                     data=nscd_all, family="poisson")
summary(meanFit1b_null)

lmtest::lrtest(meanFit1b_null, meanFit1b) # 0.00043


meanFit1d_null = lm(paste(outDiff, "~ ", controls, sep=" "),  
                    data=nscd_all)
summary(meanFit1d_null)

lmtest::lrtest(meanFit1d_null, meanFit1d) # p = 0.04545


## Adviser model
compFit1a = glm(paste(outcome, paste0("~ advHawk + presHawk", "+ formal"), sep=" "),  
                data=nscd_all, family="poisson")
summary(compFit1a)

compFit1b = glm(paste(outcome, paste0("~ advHawk + presHawk", "+", controls, "-factor(admin)")),  
                data=nscd_all, family="poisson")
summary(compFit1b)

compFit1c = lm(paste(outDiff, paste0("~ advHawk + presHawk", "+ formal"), sep=" "),  
               data=nscd_all)
summary(compFit1c)

compFit1d = lm(paste(outDiff, paste0("~ advHawk + presHawk", "+", controls, "-factor(admin)")),  
               data=nscd_all)
summary(compFit1d)


## Adviser model with administration fixed effects, without leader hawkishness
compFit2a = glm(paste(outcome, paste0("~ advHawk", "+ formal + factor(admin)"), sep=" "),  
                data=nscd_all, family="poisson")
summary(compFit2a)

compFit2b = glm(paste(outcome, paste0("~ advHawk", "+", controls)),  
                data=nscd_all, family="poisson")
summary(compFit2b)

compFit2c = lm(paste(outDiff, paste0("~ advHawk", "+ formal + factor(admin)"), sep=" "),  
               data=nscd_all)
summary(compFit2c)

compFit2d = lm(paste(outDiff, paste0("~ advHawk", "+", controls)),  
               data=nscd_all)
summary(compFit2d)


#### Table 3: Effect of Participant Hawkishness on Foreign Policy Decisions ####
stargazer(meanFit1b, meanFit1d, leaderFit1b, leaderFit1d,
          compFit1b, compFit1d, compFit2b, compFit2d,
          column.sep.width = "0pt", no.space=T, align=T, df=F,
          omit.stat = c("ll", "rsq", "adj.rsq", "aic", "f", "ser"),
          omit = c("mention", "admin"),
          covariate.labels = c("Mean Hawkishness", 
                               "Advisers' Hawkishness (Acts)", 
                               "President's Hawkishness", 
                               "No. of Attendees", 
                               "Defense", 
                               "Intelligence", 
                               "Military", 
                               "State", 
                               "Diplomatic Experience", 
                               "Intelligence Experience", 
                               "Military Experience", 
                               "5-Year MID Challenges", 
                               "US CINC",
                               "Formal"))


#### Figure 7: Summary of Three Models of Trait Aggregation ####
## (Note: This figure only shows results for full models)
plotMainHawkCoefs("")
ggsave(filename="./figures/coefPlotMain.pdf", height=4, width=7, units='in')


#### Table A9: Effect of Mean Participant Hawkishness and President's Hawkishness on Foreign Policy Decisions ####
stargazer(meanFit1a, meanFit1b, meanFit1c, meanFit1d, 
          leaderFit1a, leaderFit1b, leaderFit1c, leaderFit1d,
          column.sep.width = "0pt", no.space=T, align=T, df=F,
          omit.stat = c("ll", "rsq", "adj.rsq", "aic", "f", "ser"),
          omit = c("mention", "admin"),
          covariate.labels = c("Mean Hawkishness", 
                               "No. of Attendees", 
                               "Defense", 
                               "Intelligence", 
                               "Military", 
                               "State", 
                               "Diplomatic Experience", 
                               "Intelligence Experience", 
                               "Military Experience", 
                               "5-Year MID Challenges", 
                               "US CINC",
                               "President's Hawkishness",
                               "Formal"))

#### Table A10: Effect of Adviser Hawkishness on Foreign Policy Decisions ####
stargazer(compFit1a, compFit1b, compFit1c, compFit1d,
          compFit2a, compFit2b, compFit2c, compFit2d,
          column.sep.width = "0pt", no.space=T, align=T, df=F,
          omit.stat = c("ll", "rsq", "adj.rsq", "aic", "f", "ser"),
          omit = c("mention", "admin"),
          covariate.labels = c("Advisers' Hawkishness (Acts)", 
                               "President's Hawkishness", 
                               "No. of Attendees", 
                               "Defense", 
                               "Intelligence", 
                               "Military", 
                               "State", 
                               "Diplomatic Experience", 
                               "Intelligence Experience", 
                               "Military Experience", 
                               "5-Year MID Challenges", 
                               "US CINC",
                               "Formal"))

#### Figure A6: Summary of Three Models of Trait Aggregation ####
## (Note: This figure only shows results for all models)
plotHawkCoefs("")
ggsave(filename="./figures/coefPlotAll.pdf", height=4, width=7, units='in')


#### Table 4: Predicted Number of Conflictual Decisions Towards Adversaries ####

# Function that creates predictions for minimal and maximal values
source("./code/makePredict.R")

# Emergence model
round(makePredict(meanFit1b, "meanHawk"), 3) # 0.093, 0.608
round(makePredict(meanFit1d, "meanHawk"), 3) # 0.021, 0.346 

# Leader model
round(makePredict(leaderFit1b, "presHawk"), 3) # 0.262, 0.278
round(makePredict(leaderFit1d, "presHawk"), 3) # 0.328, 0.322 

# Adviser model, Poisson
round(makePredict(compFit1b, "advHawk"), 3) # 0.136, 0.420
round(makePredict(compFit1b, "usmidschallenge5"), 3) # 0.243, 0.248
round(makePredict(compFit1b, "cinc"), 3) # 0.203, 0.315

# Adviser model, OLS
round(makePredict(compFit1d, "advHawk"), 3) # 0.039, 0.300 
round(makePredict(compFit1d, "usmidschallenge5"), 3) # 0.173, 0.177 
round(makePredict(compFit1d, "cinc"), 3) # 0.073, 0.307


#### Section 4.1.3: Selection and Robustness ####

## Leaders vs advisers; find the intraclass correlation (ICC) ##
mixFit0 = glmer.nb(nConfAdv ~ (1|admin), data=nscd_all, family="poisson")
icc(mixFit0) # 0.179

mixFit1 = glmer.nb(nConfAdv ~ presHawk + (1|admin), data=nscd_all, family="poisson")
icc(mixFit1) # 0.177

mixFit2 = glmer(paste(outcome, "~ ", controls, "-factor(admin) + (1|admin)", sep=" "),
                data=nscd_all, family="poisson")
icc(mixFit2) # 0.037


# Save the results, which are relevant in the "meetingAnalysisApp.R" file
save.image("./data/meetingAnalysis.RData")

rm(list = ls())
